

############################################################
### Added code for factor analysis and additional models ###
############################################################
library(mirt)
library(psych)
library(corrplot)
library(car)


combo_DV <- combo %>%
  select(DV_costly, DV_interests, DV_effective, DV_strength, DV_message) %>%
  mutate(DV_costly = 100 - DV_costly)


pdf(file = "Output/corrplot.pdf")
datamatrix <- cor(combo_DV, use = "pairwise.complete.obs")
corrplot(datamatrix, method="number")
dev.off()


X <- combo_DV[, -c(1)]
Y <- combo_DV[, 1]

X <- X[complete_cases(X),]

KMO(r=cor(X))

cortest.bartlett(X)

det(cor(X))

fafitfree <- fa(X,nfactors = ncol(X), rotate = "none")
n_factors <- length(fafitfree$e.values)
scree     <- data.frame(
  Factor_n =  as.factor(1:n_factors), 
  Eigenvalue = fafitfree$e.values)
ggplot(scree, aes(x = Factor_n, y = Eigenvalue, group = 1)) + 
  geom_point() + geom_line() +
  xlab("Number of factors") +
  ylab("Initial eigenvalue") +
  labs( title = "Scree Plot", 
        subtitle = "(Based on the unreduced correlation matrix)")

pdf(file = "Output/screeplot.pdf")
parallel <- fa.parallel(X)
dev.off()


# get NAs back
X <- combo_DV[, -c(1)]


fa.none <- fa(r=combo_DV, 
              nfactors = 1, 
              # covar = FALSE, SMC = TRUE,
              fm = "pa", # type of factor analysis we want to use (“pa” is principal axis factoring)
              max.iter=100, # (50 is the default, but we have changed it to 100
              rotate= "varimax") # none rotation
print(fa.none)

fa.diagram(fa.none)

facdat <- data.frame(DV_FA = fa.none$scores)


combo <- cbind(combo, facdat)
  

model_add_full_fac <- PA1 ~ Policy + Background + badbehavior1*country1 + Inaction + Author + Rationale + EffectOrd + Prior3
# model_add_full <- DV_message ~ Policy + Background + badbehavior1 + Inaction + Author + Rationale + EffectOrd + country1


amces_add_full_fac <- cj(combo, model_add_full_fac,
                     id = ~1, alpha=0.05, estimate = "amce")

amces_add_full_fac <- amces_add_full_fac %>%
  mutate(feature = dplyr::recode(feature, "badbehavior1" = "Issue", "Inaction" = "If no action", 
                          "EffectOrd" = "Effect on target economy", "country1" = "Country",
                          "Policy" = "Proposed policy", "Prior3" = "Prior attitude"),
         level = dplyr::recode(level, "terror" = "Terror", "drugs" = "Drugs", 
                        "developnukes" = "WMD", "stopdemocracy" = "Democracy crackdown",
                        "TickingClock" = "Ticking clock", "Continue" = "Continuation",
                        "Appointee_only" = "Appointee", "Appointee_donor" = "Appointee donor",
                        "Professional_expert" = "Professional expert"))

full_a_fac <- plot(rbind(amces_add_full)) + scale_colour_viridis_d( begin = 0, end = .7) +
  labs(x="Average marginal componenet effects")+ theme(legend.position="none")


mms_add_full_fac <- cj(combo, model_add_full,
                   id = ~1, alpha=0.05, estimate = "mm")

mms_add_full_fac <- mms_add_full_fac %>%
  mutate(feature = dplyr::recode(feature, "badbehavior1" = "Issue", "Inaction" = "If no action", 
                          "EffectOrd" = "Effect on target economy", "country1" = "Country",
                          "Policy" = "Proposed policy", "Prior3" = "Prior attitude"),
         level = dplyr::recode(level, "terror" = "Terror", "drugs" = "Drugs", 
                        "developnukes" = "WMD", "stopdemocracy" = "Democracy crackdown",
                        "TickingClock" = "Ticking clock", "Continue" = "Continuation",
                        "Appointee_only" = "Appointee", "Appointee_donor" = "Appointee donor",
                        "Professional_expert" = "Professional expert"))

full_m_fac <- plot(rbind(mms_add_full_fac)) + scale_colour_viridis_d( begin = 0, end = .7) +
  labs(x="Marginal means")+ theme(legend.position="none")


library("ggpubr")
ggarrange(full_a_fac, full_m_fac, ncol = 2)

ggsave("Output/Rplot_combined_full_fac.pdf",
       width = 10, height = 6, dpi = 300, units = "in", device='pdf')

#H3 and H6 (plus others subset as well)
model_bypol_full_fac <- PA1 ~ Background + Prior3+ badbehavior1*country1  + Inaction + Author + Rationale + EffectOrd


amces_bypol_full_fac <- cj(combo, model_bypol_full_fac, 
                       id = ~1, alpha=0.05, estimate = "amce", by = ~ Policy)


amces_bypol_full_fac <- amces_bypol_full_fac %>%
  mutate(feature = dplyr::recode(feature, "badbehavior1" = "Issue", "Inaction" = "If no action", "Prior3" = "Prior attitude", 
                          "EffectOrd" = "Effect on target economy", "country1" = "Country"),
         level = dplyr::recode(level, "terror" = "Terror", "drugs" = "Drugs", 
                        "developnukes" = "WMD", "stopdemocracy" = "Democracy crackdown",
                        "TickingClock" = "Ticking clock", "Continue" = "Continuation",
                        "Appointee_only" = "Appointee", "Appointee_donor" = "Appointee donor",
                        "Professional_expert" = "Professional expert"))


plot(rbind(amces_bypol_full_fac)) + ggplot2::facet_wrap(~BY, ncol = 4L) + scale_color_jco() +
  labs(x="Average marginal component effects") + theme(legend.position="none")


ggsave("Output/Rplot_bypol_full_fac.pdf",
       width = 6.5, height = 6, dpi = 300, units = "in", device='pdf')


#H3 and H6 (plus others subset as well) - marginal means
model_bypol_full_fac <- PA1 ~ Background + Prior3 + badbehavior1*country1 + Inaction + Author + Rationale + EffectOrd

mms_bypol_full_fac <- cj(combo, model_bypol_full_fac, 
                     id = ~1, alpha=0.05, estimate = "mm", by = ~ Policy)

diff_mms_mms_bypol_full_fac <- cj(combo, model_bypol_full_fac,
                              id = ~1, alpha=0.05, estimate = "mm_diff", by = ~ Policy)


mms_bypol_full_fac <- mms_bypol_full_fac %>%
  mutate(feature = dplyr::recode(feature, "badbehavior1" = "Issue", "Inaction" = "If no action", "Prior3" = "Preexisting attitude", 
                          "EffectOrd" = "Effect on target economy", "country1" = "Country"),
         level = dplyr::recode(level, "terror" = "Terror", "drugs" = "Drugs", 
                        "developnukes" = "WMD", "stopdemocracy" = "Democracy crackdown",
                        "TickingClock" = "Ticking clock", "Continue" = "Continuation",
                        "Appointee_only" = "Appointee", "Appointee_donor" = "Appointee donor",
                        "Professional_expert" = "Professional expert"))

diff_mms_mms_bypol_full_fac <- diff_mms_mms_bypol_full_fac %>%
  mutate(feature = dplyr::recode(feature, "badbehavior1" = "Issue", "Inaction" = "If no action", "Prior3" = "Preexisting attitude", 
                          "EffectOrd" = "Effect on target economy", "country1" = "Country"),
         level = dplyr::recode(level, "terror" = "Terror", "drugs" = "Drugs", 
                        "developnukes" = "WMD", "stopdemocracy" = "Democracy crackdown",
                        "TickingClock" = "Ticking clock", "Continue" = "Continuation",
                        "Appointee_only" = "Appointee", "Appointee_donor" = "Appointee donor",
                        "Professional_expert" = "Professional expert"))

plot(rbind(mms_bypol_full_fac, diff_mms_mms_bypol_full_fac)) + ggplot2::facet_wrap(~BY, ncol = 4L) + scale_color_jco() +
  labs(x="Marginal means") + theme(legend.position="none")

ggsave("Output/Rplot_mm_bypol_full_fac.pdf",
       width = 6.5, height = 6, dpi = 300, units = "in", device='pdf')

    